require(mvtnorm)
require(OpenMx)
require(psych)
require(MASS)

inter <- 0                                                  # Generate the interaction parameter. We used 0, .5, 1, 1.5 & 2
nsim <- 1000                                              # Set the number of simulations (We used 100000)

AllMZdata <- list(NULL)
AllDZdata <- list(NULL)

MZdatasets <- list(NULL)
DZdatasets <- list(NULL)


nv <- 1                                                     # Set the number of variables that you will be using (per person)
Vars <- c('t')                                              # Name the variable
selVars <- paste(Vars,c(rep(1,nv),rep(2,nv)),sep="")        # Create an object that will select the variables
ntv <- nv*2                                                 # Set the total number of variable that you will use (1 var * 2 twins)


mzmat <- matrix(c(1,1,0,0,0,0,
                  1,1,0,0,0,0,
                  0,0,1,1,0,0,
                  0,0,1,1,0,0,
                  0,0,0,0,1,0,
                  0,0,0,0,0,1),6)                           # Create a matrix of correlations for the genetic and environmental variance components for each twin (MZ)

dzmat <- matrix(c(1,1/2,0,0,0,0,
                  1/2,1,0,0,0,0,
                  0,0,1,1,0,0,
                  0,0,1,1,0,0,
                  0,0,0,0,1,0,
                  0,0,0,0,0,1),6)                           # Create a matrix of correlations for the genetic and environmental variance components for each twin (DZ)



Avar <- (0:6/10)                                            # Create an object with each of the possible Additive Genetic Parameters
Cvar <- (6:0/10)                                            # Create an object with each of the possible Common Environmental Parameters


for(param in 1:7){
 for(sim in 1:nsim) {

	 mzSim  <- mvrnorm(n=1000, c(0,0,0,0,0,0), mzmat , empirical=T)   # Simulate genetic and environmental factor scores for each MZ twin
	 dzSim  <- mvrnorm(n=1000, c(0,0,0,0,0,0), dzmat , empirical=T)   # Simulate genetic and environmental factor scores for each DZ twin

	 VarCompA <- Avar[param]                                          # Set the Genetic Variance Component
	 VarCompC <- Cvar[param]                                          # Set the Common Environmental Variance Component
	 VarCompE <- .4                                                   # Set the Unique Environmental Variance Component

	 aPath <- sqrt(VarCompA)                                          # Set the Genetic Variance Path Coefficeint
	 cPath <- sqrt(VarCompC)                                          # Set the Common Environmental Path Coefficeint
	 ePath <- sqrt(VarCompE)                                          # Set the Unique Environmental Path Coefficeint



	 MZdata <- cbind((mzSim[,1]*aPath+
	                  mzSim[,3]*cPath+
	                  mzSim[,5]*ePath+
	                  mzSim[,1]*mzSim[,3]*inter),
	                 (mzSim[,2]*aPath+
	                  mzSim[,4]*cPath+
	                  mzSim[,6]*ePath+
	                  mzSim[,2]*mzSim[,4]*inter))                      # Construct the phenotype by adding the genetic/environmental factor score * the path + the bias due to the interaction for MZ twins

	 DZdata <- cbind((dzSim[,1]*aPath+
	                  dzSim[,3]*cPath+
	                  dzSim[,5]*ePath+
	                  dzSim[,1]*dzSim[,3]*inter),
	                 (dzSim[,2]*aPath+
	                  dzSim[,4]*cPath+
	                  dzSim[,6]*ePath+
	                  dzSim[,2]*dzSim[,4]*inter))                      # Construct the phenotype by adding the genetic/environmental factor score * the path + the bias due to the interaction for DZ twins

	 colnames (MZdata) <- selVars                                       # Give the datasets column names
	 colnames (DZdata) <- selVars

		
		AllMZdata[[sim]] <- MZdata	
		AllDZdata[[sim]] <- DZdata	
		
	}
	MZdatasets[[param]]<- AllMZdata
	DZdatasets[[param]]<- AllDZdata
}


#MZdatasets
#DZdatasets
